﻿/* Copyright 2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Marcus Cuda
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    /// <summary>
    /// Computes the QR decomposition of a <see cref="Matrix"/> using the Householder method.
    /// </summary>
    internal abstract class AbstractHouseholder
    {
        protected readonly int mColumns;
        protected readonly int mRows;
        private bool mComputed;
        protected double mDeterminant = 1;
        protected bool mIsFullRank = true;
        protected Matrix mQ;
        protected Matrix mR;


        /// <summary>
        /// Constructs the QR object for the given matrix.
        /// </summary>
        /// <param name="matrix">The matrix to decompose.</param>
        protected AbstractHouseholder(Matrix matrix)
        {
            mR = matrix.Clone();
            mRows = matrix.Rows;
            mColumns = matrix.Columns;
        }

        /// <summary>
        /// Returns the orthogonal Q matrix.
        /// </summary>
        /// <returns>The orthogonal Q matrix.</returns>
        public Matrix Q()
        {
            Compute();
            return mQ.Clone();
        }

        /// <summary>
        /// Returns the upper triangular factor R.
        /// </summary>
        /// <returns>The upper triangular factor R.</returns>
        public Matrix R()
        {
            Compute();
            return mR.GetUpperTriangle();
        }

        /// <summary>
        /// Determines whether the matrix is full rank or not.
        /// </summary>
        /// <returns><b>true</b> if the matrix is full rank; otherwise <b>false</b>.</returns>         
        public bool IsFullRank()
        {
            Compute();
            return mIsFullRank;
        }

        ///<summary>Calculates the determinant (absolute value) of the matrix.</summary>
        ///<returns>The determinant of the matrix.</returns>
        ///<exception cref="MatrixNotSquareException">If the matrix is not square.</exception>
        public double Determinant()
        {
            Compute();
            return mDeterminant;
        }

        private void Compute()
        {
            if (!mComputed)
            {
                DoCompute();
                for (int i = 0; i < mR.Columns; i++)
                {
                    mDeterminant *= mR.ValueAt(i, i);
                    if (Precision.EqualsWithinDecimalPlaces(System.Math.Abs(mR.ValueAt(i, i)), 0.0, Constants.DecimalsForComparison))
                    {
                        mIsFullRank = false;
                        mDeterminant = 0;
                        break;
                    }
                }
                mDeterminant = System.Math.Abs(mDeterminant);
                mComputed = true;
            }
        }

        protected abstract void DoCompute();

        public Matrix Solve(Matrix input)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (mRows != input.Rows)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            Compute();
            Matrix result = mR.CreateMatrix(mR.Columns, input.Columns);
            Solve(input, result);
            return result;
        }

        public void Solve(Matrix input, Matrix result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (input.Columns != result.Columns)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (mRows != input.Rows)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mColumns != result.Rows)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }
            Compute();

            Matrix b = input.Clone();
            DoSolve(b, result);
        }

        public Vector Solve(Vector input)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (mRows != input.Count)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            Compute();
            Vector result = mQ.CreateVector(mR.Columns);
            Solve(input, result);
            return result;
        }

        public void Solve(Vector input, Vector result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (mRows != input.Count)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mColumns != result.Count)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            Compute();

            Vector b = input.Clone();
            DoSolve(b, result);
        }

        protected abstract void DoSolve(Matrix b, Matrix result);

        protected abstract void DoSolve(Vector b, Vector result);

        protected static double Norm(double[] vector)
        {
            double ret = 0.0;
            for (int i = 0; i < vector.Length; ++i)
            {
                ret += vector[i]*vector[i];
            }
            return System.Math.Sqrt(ret);
        }
    }
}
